---
title: "Insecurity and political values in the Arab world"
author: "Melani Cammett, Ishac Diwan & Irina Vartanova"
output: html_document
---

```{r setup, include = FALSE}

knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
                      dpi = 300)

library(tidyverse)
library(cowplot)
library(lfe)
library(broom)
library(stargazer)

wvs <- read_csv("wvs-aggr.csv")
ab <- read_csv("ab-security-and-values.csv",
               col_types = cols(
                 .default = col_double(),
                 country = col_character(),
                 sex = col_factor(),
                 age_gr = col_factor(),
                 edu = col_factor(),
                 income = col_factor()))


ab <- ab %>% 
  mutate(cntry_year = str_c(country, "_", year),
          year_fct = as.factor(year),
          age_gr = factor(age_gr, 
                         levels = c("18-29", "30-49", "50>"))) %>% 
  group_by(cntry_year) %>% 
  mutate(econ_sec_natl = mean(econ_sec, na.rm = TRUE),
         pers_sec_natl = mean(pers_sec, na.rm = TRUE),
         econ_sec_ind = econ_sec - econ_sec_natl, 
         pers_sec_ind = pers_sec - pers_sec_natl,
         islam_natl = mean(islamism , na.rm = TRUE),
         islam_ind = islamism  - islam_natl) 

selected_cntries <- c("Egypt", "Tunisia", "Morocco", "Jordan", "Algeria", "Iraq")

ab6 <- filter(ab, country %in% selected_cntries)

```

## Figures

### Figure 1

```{r, fig.width = 7, fig.asp = .8, fig.cap= "Figure 1. (a–d) Economic and personal insecurity, preference for democracy, and trust in government, in the Arab countries and around the world. Source: World Values Survey, Wave 6; for GDP per capita, World Development Indicators 2015."}

pl_aggr <- wvs %>% 
  gather("variable", "value", econ_sec, pers_sec, trust_gov, pfd) %>% 
  mutate(variable = factor(variable, 
                           levels = c("econ_sec", "pers_sec",
                                       "trust_gov", "pfd"),
                           labels = c("Economic security", 
                                      "Personal security",
                                      "Trust in Government",
                                      "Commitment to Democracy"
                                      )))



pl <- pl_aggr %>% 
  ggplot(aes(log(gdpc), value, label = cntry_code, color = arab))+
  geom_point(size = 2) + 
  ggrepel::geom_text_repel(size = 2) +
  scale_x_continuous()+
  facet_wrap(~ variable, scales = "free_y")+
  scale_color_manual(values = c("grey70", "grey10")) + 
  labs(x = "log GDPc", y = NULL, color = NULL)+
  theme_cowplot() +
  theme(legend.position = c(.85, .6))


pl <- ggdraw(pl) +
  draw_plot_label(c("a", "b", "c", "d"), 
            size = 16,
            x = c(.015, .51, .015, .51), 
            y = c(.995, .995, .545, .545))
pl

# save_plot("fig1_security-and-values-by-gdp.eps", pl, base_width = 7, base_height = 5.6)

```


### Figure 2


```{r, fig.width = 14, fig.asp = .3, fig.cap = "Figure 2. Cross-time values of security and political values in selected Arab countries. Source: Arab Barometer, Waves 1–4."}

trends_data <- ab6 %>% 
  group_by(country, year) %>% 
  summarise(econ_sec = Hmisc::wtd.mean(econ_sec, wt),
            pers_sec = Hmisc::wtd.mean(pers_sec, wt),
            demo_pref = Hmisc::wtd.mean(demo_pref, wt),
            trust_governm = Hmisc::wtd.mean(trust_governm, wt)) %>% 
  gather(variable, value, econ_sec:trust_governm) %>% 
  mutate(variable = factor(variable, 
                           levels = c("econ_sec", "pers_sec", 
                                      "trust_governm", "demo_pref"),
                           labels = c("Economic security", 
                                      "Personal security",
                                      "Trust in Govrnment",
                                      "Commitment to Democracy")),
         label = sprintf("%.2f", value)) 

text_left  <- trends_data %>% 
  group_by(variable, country) %>% 
  filter(year == min(year)) %>% 
  mutate(label = sprintf("%.2f", value))

text_right  <- trends_data %>% 
  group_by(variable, country) %>% 
  filter(year == max(year)) %>% 
  mutate(label = sprintf("%.2f", value))

pl2 <- trends_data %>% 
  ggplot(aes(year, value, label = label)) +
  geom_line(color = "midnightblue") +
  geom_text(data = text_left, color = "midnightblue", nudge_x = -1.2, nudge_y = .01, size = 2.5) +
  geom_text(data = text_right, color = "midnightblue", nudge_x = 1.2, nudge_y = .01, size = 2.5) +
  facet_grid(variable ~ country, switch = "y")+ #, scales = "free_y")+
  scale_x_continuous(limits = c(2005, 2018), 
                     breaks = c(2006, 2011, 2016))+
  scale_y_continuous(position = "right")+
  labs(x = NULL, y = NULL)+
  theme_cowplot() +
  theme(axis.text.y = element_text(size = 11),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 10),
        strip.background.y = element_blank(),
        strip.text.y = element_text(angle = 180, hjust = 1, margin = margin(0, 10, 0, 0))) 

pl2
# save_plot("fig2_cross-time-values.eps", pl2, base_width = 10, base_height = 3.5)

```


### Figure 3

```{r, fig.width = 7, fig.asp = .5, fig.cap="Figure 3. (a and b) Evolution of economic insecurity and political values in selected Arab countries over time. Source: Arab Barometer, Waves 1–4. Source: Arab Barometer, Waves 1–4."}

pl3 <- ab6 %>% 
  group_by(country, year) %>% 
  summarise(econ_sec_natl = Hmisc::wtd.mean(econ_sec, wt),
            demo_pref_natl = Hmisc::wtd.mean(demo_pref, wt),
            trust_governm_natl = Hmisc::wtd.mean(trust_governm, wt)) %>% 
  # add WVS for info on country codes
  left_join(wvs %>% select(country = cntry, cntry_code)) %>% 
  gather(variable, value, demo_pref_natl, trust_governm_natl) %>% 
  mutate(cntry_year = str_c(cntry_code, " ", 
                            str_replace(year, "^20", "")),
         variable = factor(variable, 
                           labels = c("Commitment to Democracy",
                                      "Trust in Government"))) %>% 
  ggplot(aes(econ_sec_natl, value, label = cntry_year)) +
  geom_point(color = "midnightblue") +
  ggrepel::geom_text_repel(size = 3, color = "grey30") +
  theme_cowplot() + 
  labs(x = "Economic security", y = NULL, color = NULL) +
  geom_smooth(method = "lm", se = FALSE, color = "grey50") +
  facet_wrap(~variable, scales = "free")

pl3 <- ggdraw(pl3) +
  draw_plot_label(c("a", "b"), 
            size = 16,
            x = c(.015, .51), 
            y = c(.99, .99))
pl3

# save_plot("fig3_evolution-selected-cntries.eps", pl3, base_width = 7, base_height = 3.5)


```



## Tables

### Table 1

```{r descriptive}

cont <- c("demo_pref", "trust_governm", 
         "econ_sec", "econ_sec_natl", 
         "pers_sec", "pers_sec_natl",
         "islam_ind ", "islam_year")

factors <- c("sex", "age_gr", "edu", "income")
order <- c(levels(ab$sex), 
           levels(ab$age_gr), 
           levels(ab$edu), 
           levels(ab$income))

repl_pattern <- c("demo_pref" = "PD", 
                  "trust_governm" = "TD", 
                  "econ_sec$" = "ES individual",
                  "econ_sec_natl" = "ES country/time", 
                  "pers_sec$" = "PS individual", 
                  "pers_sec_natl" = "PS country/time",
                  "islam_year" = "PI country/time", 
                  "islam_ind " = "PI individual",
                  "sex" = "",
                  "age_gr" = "Age ",
                  "edu" = "Edu-",
                  "income" = "Income-")

descr_13 <- ab %>% 
  ungroup() %>% 
  select(one_of(cont), wt) %>% 
  gather(variable, value, -wt) %>% 
  group_by(variable) %>% 
  summarise(mean = Hmisc::wtd.mean(value, wt),
            sd = sqrt(Hmisc::wtd.var(value, wt)), 
            min = min(value, na.rm = TRUE),
            max = max(value, na.rm = TRUE)) %>% 
  mutate(variable = factor(variable, levels = cont)) %>% 
  arrange(variable) %>% 
  bind_rows(ab %>% 
              ungroup() %>% 
              select(one_of(factors), wt) %>% 
              gather(variable, value, -wt) %>% 
              drop_na() %>% 
              group_by(variable) %>% 
              count(value, wt = wt) %>% 
              mutate(mean = n/sum(n)) %>% 
              ungroup() %>% 
              mutate(value = factor(value, levels = order),
                     variable = str_c(variable, value, " ")) %>% 
              arrange(value) %>% 
              select(-n, -value)) %>% 
  mutate(variable = str_replace_all(variable, repl_pattern)) %>% 
  mutate_at(vars(-variable), ~sprintf("%.2f", .) %>% 
              str_replace("NA", "-")) %>%
  rename_at(vars(-variable), ~str_c(., "_13"))
  
descr_6 <- ab6 %>% 
  ungroup() %>% 
  select(one_of(cont), wt) %>% 
  gather(variable, value, -wt) %>% 
  group_by(variable) %>% 
  summarise(mean = Hmisc::wtd.mean(value, wt),
            sd = sqrt(Hmisc::wtd.var(value, wt)), 
            min = min(value, na.rm = TRUE),
            max = max(value, na.rm = TRUE)) %>% 
  mutate(variable = factor(variable, levels = cont)) %>% 
  arrange(variable) %>% 
  bind_rows(ab6 %>% 
              ungroup() %>% 
              select(one_of(factors), wt) %>% 
              gather(variable, value, -wt) %>% 
              drop_na() %>% 
              group_by(variable) %>% 
              count(value, wt = wt) %>% 
              mutate(mean = n/sum(n)) %>% 
              ungroup() %>% 
              mutate(value = factor(value, levels = order),
                     variable = str_c(variable, value, " ")) %>% 
              arrange(value) %>% 
              select(-n, -value)) %>% 
  mutate(variable = str_replace_all(variable, repl_pattern)) %>% 
  mutate_at(vars(-variable), ~sprintf("%.2f", .) %>% 
              str_replace("NA", "-")) %>%
  rename_at(vars(-variable), ~str_c(., "_6"))

descr_6 %>% 
  left_join(descr_13) %>% 
  knitr::kable(col.names = c("", rep(c("Mean or %", "SD", "Min", "Max"), 2))) %>% 
  kableExtra::add_header_above(c(" " = 1, "6 countrie" = 4, "13 countries" = 4))

```



### Table 2

```{r, results="asis"}

dp1 <- felm(demo_pref ~ econ_sec_ind + econ_sec_natl  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

dp2 <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl  + islam_natl + islam_ind  + 
             sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

tg1 <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

tg2 <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

stargazer(dp1, dp2, tg1, tg2, 
          type = "html", 
          report = "vc*s",
          star.cutoffs = c(.05, .01, .001),
          no.space = TRUE,
          single.row = FALSE,
          column.labels = c("Preference for democracy",
                            "Trust in Government"),
          column.separate = c(2, 2),
          title = "Table A2. Perceptions of insecurity by population subgroups in the Arab region, 6 selected countries.")

```

## Supplemental online appendix

### Table A2

```{r, results = "asis"}

ec_sec1 <- felm(econ_sec ~ islam_natl + islam_ind  + sex + age_gr + edu + 
                    income | country | 0 | cntry_year,
                ab6, 
                weights = ab6$wt)

ec_sec2 <- felm(econ_sec ~  pers_sec + islam_natl + islam_ind  + sex + age_gr + edu + 
                  income | country | 0 | cntry_year,
                ab6, 
                weights = ab6$wt)

pers_sec1 <- felm(pers_sec ~  islam_natl + islam_ind  + sex + age_gr + edu + 
                      income| country | 0 | cntry_year,
                ab6, 
                weights = ab6$wt)

pers_sec2 <- felm(pers_sec ~  econ_sec + islam_natl + islam_ind  +  sex + age_gr + edu + 
                  income| country | 0 | cntry_year,
                ab6, 
                weights = ab6$wt)

stargazer::stargazer(ec_sec1, ec_sec2, pers_sec1, pers_sec2, 
                     type = "html", 
                     column.labels = c("Economic Security",
                                       "Personal Security"),
                     column.separate = c(2, 2),
                     title = "Table A2. Perceptions of insecurity by population subgroups in the Arab region, 6 selected countries.",
                     report = "vc*s",
                     star.cutoffs = c(.05, .01, .001),
                     no.space = TRUE,
                     single.row = FALSE)
```



### Table A3

```{r, results = "asis"}

ec_sec1_12 <- felm(econ_sec ~  islam_natl + islam_ind  + sex + age_gr + edu + 
                     income | country | 0 | cntry_year,
                ab, 
                weights = ab$wt)

ec_sec2_12 <- felm(econ_sec ~ pers_sec + islam_natl + islam_ind  + sex + age_gr + edu + 
                  income  | country | 0 | cntry_year,
                ab, 
                weights = ab$wt)

pers_sec1_12 <- felm(pers_sec ~  islam_natl + islam_ind  + sex + age_gr + edu + 
                       income  | country | 0 | cntry_year,
                ab, 
                weights = ab$wt)

pers_sec2_12 <- felm(pers_sec ~  econ_sec + islam_natl + islam_ind  + sex + age_gr + edu + 
                  income  | country | 0 | cntry_year,
                ab, 
                weights = ab$wt)

stargazer::stargazer(ec_sec1_12, ec_sec2_12, 
                     pers_sec1_12, pers_sec2_12,
                     type = "html", 
                     column.labels = c("Economic Security",
                                       "Personal Security"),
                     column.separate = c(2, 2),
                     title = "Table A3. Perceptions of insecurity, thirteen countries.",
                     report = "vc*s",
                     star.cutoffs = c(.05, .01, .001),
                     no.space = TRUE,
                     single.row = FALSE)
```

### Table A4

```{r, results = "asis"}

dp1_12 <- felm(demo_pref~ econ_sec_ind + econ_sec_natl + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab, 
           weights = ab$wt)

dp2_12 <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab, 
           weights = ab$wt)


tg1_12 <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab, 
           weights = ab$wt)

tg2_12 <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab, 
           weights = ab$wt)

stargazer::stargazer(dp1_12, dp2_12, 
                     tg1_12, tg2_12, 
                     type = "html", 
                     column.labels = c("Economic Security",
                                       "Personal Security"),
                     column.separate = c(2, 2),
                     title = "Table A4. Perceptions of insecurity, thirteen countries.",
                     report = "vc*s",
                     star.cutoffs = c(.05, .01, .001),
                     no.space = TRUE,
                     single.row = FALSE)
```


### Table A5

```{r, results = "asis"}

econ_sec_mod <- function(variable) {
  my_formula <- str_c(variable, " ~ econ_sec_ind + econ_sec_natl  + islam_natl + islam_ind  +
            sex + age_gr + edu + income | country | 0 | cntry_year")
  
  felm(formula = as.formula(my_formula), 
       ab6, 
       weights = ab6$wt)
}

pers_sec_mod <- function(variable) {
  my_formula <- str_c(variable, " ~ pers_sec_ind + pers_sec_natl + islam_natl + islam_ind  +
            sex + age_gr + edu + income | country | 0 | cntry_year")
  
  felm(formula = as.formula(my_formula), 
       ab6, 
       weights = ab6$wt)
}

trust_vars <- str_subset(names(ab6), "trust_")
trust_vars <- str_subset(trust_vars, "trust_(governm|mb)", negate = TRUE)

trust_models <- c(map(trust_vars, econ_sec_mod), map(trust_vars, pers_sec_mod))


stargazer::stargazer(trust_models, 
                     type = "html", 
                     report = "vc*s",
                     column.labels = c("courts", "parlm",  
                                       "police", "army",   
                                       "courts", "parlm",  
                                       "police", "army"),
                     title = "Table A5. Disaggregated Trust in Government variables, six countries.",
                    star.cutoffs = c(.05, .01, .001),
                     no.space = TRUE,
                     single.row = FALSE)
```


### Table A6

```{r results = "asis"}

dp_vars <- c(str_subset(names(ab6), "q516"), 
              "demo_auto", "democr_approval")

dp_models <- c(map(dp_vars, econ_sec_mod), 
                map(dp_vars, pers_sec_mod))


stargazer::stargazer(dp_models,
                     type = "html", 
                     title = "Table A6. Alternative measures of Preference for Democracy, six countries.",
                     column.labels = c("q5161", "q5162", "q5163",
                                       "Demo-auto", "Demo approval",
                                       "q5161", "q5162", "q5163",
                                       "Demo-auto", "Demo approval"),
                     report = "vc*s",
                     star.cutoffs = c(.05, .01, .001),
                     no.space = TRUE,
                     single.row = FALSE)
```


### Table A7 

```{r, results = "asis"}

by_cntry <- ab6 %>% 
  ungroup() %>% 
  mutate(country = as.character(country)) %>% 
  group_by(country) %>% 
  nest()
  

dp_by_cntry <- by_cntry %>% 
  mutate(dp_ec = map(data,
                      ~lm(demo_pref ~ econ_sec + islamism  + sex + age_gr + edu + 
                            income + year_fct,
                          data = ., weights = wt)),
         dp_pers = map(data,
                        ~lm(demo_pref ~ pers_sec + islamism  + sex + age_gr + edu + 
                              income + year_fct,
                            data = ., weights = wt))) %>% 
  gather(variable, model, dp_ec:dp_pers) %>% 
  arrange(country)

stargazer(dp_by_cntry$model, 
          type = "html", 
          title = "Table A7. Country level regressions: Preference for Democracy.",
          column.labels = dp_by_cntry$country,
          report = "vc*s",
          star.cutoffs = c(.05, .01, .001),
          no.space = TRUE,
          single.row = FALSE)

```

### Table A8

```{r, results = "asis"}

tg_by_cntry <- by_cntry %>% 
  mutate(tg_ec = map(data,
                     ~lm(trust_governm ~ econ_sec + islamism  + sex + age_gr + edu + 
                           income + year_fct,
                         data = ., weights = wt)),
         tg_pers = map(data,
                       ~lm(trust_governm ~ pers_sec + islamism  + sex + age_gr + edu + 
                             income + year_fct,
                           data = ., weights = wt))) %>% 
   gather(variable, model, tg_ec:tg_pers) %>% 
  arrange(country)

stargazer(tg_by_cntry$model, 
          type = "html", 
          title = "Table A8. Country level regressions: Trust in government.",
          column.labels = tg_by_cntry$country,
          report = "vc*s",
          star.cutoffs = c(.05, .01, .001),
          no.space = TRUE,
          single.row = FALSE)

```

### Table A9. The Effects of Security on Political Values by Population Subgroups

```{r values_models}


dp_econ_edu <- felm(demo_pref ~ econ_sec_ind + econ_sec_natl*edu  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

dp_econ_inc <- felm(demo_pref ~ econ_sec_ind + econ_sec_natl*income  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
dp_econ_age <- felm(demo_pref ~ econ_sec_ind + econ_sec_natl*age_gr  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
dp_econ_isl <- felm(demo_pref ~ econ_sec_ind + econ_sec_natl*islam_ind   + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)


dp_pers_edu <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl*edu  + islam_natl + islam_ind  + 
             sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

dp_pers_inc <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl*income  + islam_natl + islam_ind  + 
             sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

dp_pers_age <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl*age_gr  + islam_natl + islam_ind  + 
             sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

dp_pers_isl <- felm(demo_pref ~ pers_sec_ind + pers_sec_natl*islam_ind   + islam_natl + islam_ind  + 
             sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)


tg_econ_edu <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl*edu  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
tg_econ_inc <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl*income  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
tg_econ_age <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl*age_gr  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

tg_econ_isl <- felm(trust_governm ~ econ_sec_ind + econ_sec_natl*islam_ind   + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

tg_pers_edu <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl*edu  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
tg_pers_inc <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl*income  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)
tg_pers_age <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl*age_gr  + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)

tg_pers_isl <- felm(trust_governm ~ pers_sec_ind + pers_sec_natl*islam_ind   + islam_natl + islam_ind  + 
            sex + age_gr + edu + income | country | 0 | cntry_year, 
           ab6, 
           weights = ab6$wt)


format_mod <- function(mod, int_term){
  tidy(mod) %>% 
    filter(str_detect(term, int_term)) %>% 
    mutate(stars = case_when(
      p.value >.05 ~ "",
      p.value < .001 ~ "***",
      p.value < .01 ~ "**",
      p.value < .05 ~ "*"
    )) %>% 
    transmute(term, 
              coef = sprintf("%.3f%s ", estimate, stars),
              se = sprintf("(%.3f)", std.error)) %>% 
    gather(type, coef, coef, se) %>% 
    mutate(term = factor(term, levels = unique(term))) %>% 
    arrange(term)
}

dp_econ_int <- bind_rows(
  format_mod(dp_econ_edu, "edu|sec_natl"), 
  format_mod(dp_econ_age, "age|sec_natl"),
  format_mod(dp_econ_inc, "inc|sec_natl"),
  format_mod(dp_econ_isl, "islam_ind|sec_natl")
  ) 

tg_econ_int <- bind_rows(
  format_mod(tg_econ_edu, "edu|sec_natl"),
  format_mod(tg_econ_age, "age|sec_natl"),
  format_mod(tg_econ_inc, "inc|sec_natl"),
  format_mod(tg_econ_isl, "islam_ind|sec_natl")
  ) 

dp_pers_int <- bind_rows(
  format_mod(dp_pers_edu, "edu|sec_natl"),
  format_mod(dp_pers_age, "age|sec_natl"),
  format_mod(dp_pers_inc, "inc|sec_natl"),
  format_mod(dp_pers_isl, "islam_ind|sec_natl")
) 

tg_pers_int <-   bind_rows(
  format_mod(tg_pers_edu, "edu|sec_natl"),
  format_mod(tg_pers_age, "age|sec_natl"),
  format_mod(tg_pers_inc, "inc|sec_natl"),
  format_mod(tg_pers_isl, "islam_ind|sec_natl")
) 

bind_cols(
  dp_econ_int %>% select(term1 = term, dp_econ = coef),
  tg_econ_int %>% select(tg_econ = coef),
  dp_pers_int %>% select(term2 = term, dp_pers = coef),
  tg_pers_int %>% select(tg_pers = coef)
) %>% 
  knitr::kable() %>% 
  kableExtra::collapse_rows()

```

